Analyse tree
David Helekal
12/2/2020
Real-life phylogeny analysis
In this notebook a single tree read from a newick string is analysed. We proceed by inferring the posterior via RjMCMC, followed by it’s analyses.
MCMC Controls
First we set up the MCMC parameters
set.seed(2)
burn_in <- 0.2
n_it <- 1e8
thinning <- n_it/1e4Priors
Next, we parametrise and ready the priors for inference. We use the following parameters: - \(poisson(1)\) as the prior on the number of expansions - \(lognormal(3,3)\) as a relatively uninformative prior on parent population size - \(exp(10)\) as the prior on how long it takes for an expansion to reach half its carrying capacity in non-dimensional time units - \(lognormal(\log(N), 1/2)\) where \(N\) is the parent population size, as the prior on carrying capacities of individual expansions - \(gamma\) with mean \(\frac{N}{2}\) and variance \(\frac{N}{2}\) where \(N\) is the parent population size as the prior on when expansions happen N.B.: One should note that the model is conditional on the size of the parent population which effectively determines the tree scale. As such it is imperative that there is sufficient signal to infer this parameter with high enough confidence.
priors <- standard_priors(expansion_rate=1,
N_mean_log=3,
N_sd_log=3,
t_mid_rate=5,
K_sd_log=1,
exp_time_nu=1/2,
exp_time_kappa=1/2)Run MCMC
We now read the tree and run the MCMC.
tree_in <- list.files(pattern=".*\\.nwk")
if (length(tree_in) > 0) {
tree <- read.tree(file = tree_in)
} else {
tree_in <- list.files(pattern=".*\\.tre")
tree <- read.tree(file = tree_in)
}
tree <- makeNodeLabel(tree)
tree <- ladderize(tree, right = F)
if(length(tree$edge.length[tree$edge.length<0]) > 0) print("Negative branch lengths encountered and rounded to 1e-8")
tree$edge.length[tree$edge.length<0]<-1e-8
start <- proc.time()
expansions <- run_expansion_inference(tree, priors, 1, n_it=n_it, thinning=thinning)##
|
| | 0%
|
| | 1%
|
|= | 1%[1] "The acceptance rate is: 0.2108"
##
|
|= | 2%[1] "The acceptance rate is: 0.2083"
##
|
|== | 2%
|
|== | 3%[1] "The acceptance rate is: 0.2067"
##
|
|== | 4%
|
|=== | 4%[1] "The acceptance rate is: 0.2059"
##
|
|=== | 5%
|
|==== | 5%[1] "The acceptance rate is: 0.2052"
##
|
|==== | 6%[1] "The acceptance rate is: 0.2052"
##
|
|===== | 6%
|
|===== | 7%[1] "The acceptance rate is: 0.2051"
##
|
|===== | 8%
|
|====== | 8%[1] "The acceptance rate is: 0.2048"
##
|
|====== | 9%[1] "The acceptance rate is: 0.2046"
##
|
|======= | 9%
|
|======= | 10%[1] "The acceptance rate is: 0.205"
##
|
|======= | 11%
|
|======== | 11%[1] "The acceptance rate is: 0.2052"
##
|
|======== | 12%[1] "The acceptance rate is: 0.2054"
##
|
|========= | 12%
|
|========= | 13%[1] "The acceptance rate is: 0.2058"
##
|
|========= | 14%
|
|========== | 14%[1] "The acceptance rate is: 0.206"
##
|
|========== | 15%[1] "The acceptance rate is: 0.2064"
##
|
|=========== | 15%
|
|=========== | 16%[1] "The acceptance rate is: 0.2065"
##
|
|============ | 16%
|
|============ | 17%[1] "The acceptance rate is: 0.2064"
##
|
|============ | 18%
|
|============= | 18%[1] "The acceptance rate is: 0.2065"
##
|
|============= | 19%[1] "The acceptance rate is: 0.2065"
##
|
|============== | 19%
|
|============== | 20%[1] "The acceptance rate is: 0.2065"
##
|
|============== | 21%
|
|=============== | 21%[1] "The acceptance rate is: 0.2065"
##
|
|=============== | 22%[1] "The acceptance rate is: 0.2065"
##
|
|================ | 22%
|
|================ | 23%[1] "The acceptance rate is: 0.2065"
##
|
|================ | 24%
|
|================= | 24%[1] "The acceptance rate is: 0.2064"
##
|
|================= | 25%
|
|================== | 25%[1] "The acceptance rate is: 0.2063"
##
|
|================== | 26%[1] "The acceptance rate is: 0.2063"
##
|
|=================== | 26%
|
|=================== | 27%[1] "The acceptance rate is: 0.2064"
##
|
|=================== | 28%
|
|==================== | 28%[1] "The acceptance rate is: 0.2065"
##
|
|==================== | 29%[1] "The acceptance rate is: 0.2065"
##
|
|===================== | 29%
|
|===================== | 30%[1] "The acceptance rate is: 0.2066"
##
|
|===================== | 31%
|
|====================== | 31%[1] "The acceptance rate is: 0.2066"
##
|
|====================== | 32%[1] "The acceptance rate is: 0.2067"
##
|
|======================= | 32%
|
|======================= | 33%[1] "The acceptance rate is: 0.2066"
##
|
|======================= | 34%
|
|======================== | 34%[1] "The acceptance rate is: 0.2065"
##
|
|======================== | 35%[1] "The acceptance rate is: 0.2064"
##
|
|========================= | 35%
|
|========================= | 36%[1] "The acceptance rate is: 0.2063"
##
|
|========================== | 36%
|
|========================== | 37%[1] "The acceptance rate is: 0.2063"
##
|
|========================== | 38%
|
|=========================== | 38%[1] "The acceptance rate is: 0.2063"
##
|
|=========================== | 39%[1] "The acceptance rate is: 0.2064"
##
|
|============================ | 39%
|
|============================ | 40%[1] "The acceptance rate is: 0.2065"
##
|
|============================ | 41%
|
|============================= | 41%[1] "The acceptance rate is: 0.2065"
##
|
|============================= | 42%[1] "The acceptance rate is: 0.2066"
##
|
|============================== | 42%
|
|============================== | 43%[1] "The acceptance rate is: 0.2066"
##
|
|============================== | 44%
|
|=============================== | 44%[1] "The acceptance rate is: 0.2066"
##
|
|=============================== | 45%
|
|================================ | 45%[1] "The acceptance rate is: 0.2066"
##
|
|================================ | 46%[1] "The acceptance rate is: 0.2067"
##
|
|================================= | 46%
|
|================================= | 47%[1] "The acceptance rate is: 0.2066"
##
|
|================================= | 48%
|
|================================== | 48%[1] "The acceptance rate is: 0.2066"
##
|
|================================== | 49%[1] "The acceptance rate is: 0.2065"
##
|
|=================================== | 49%
|
|=================================== | 50%[1] "The acceptance rate is: 0.2065"
##
|
|=================================== | 51%
|
|==================================== | 51%[1] "The acceptance rate is: 0.2064"
##
|
|==================================== | 52%[1] "The acceptance rate is: 0.2063"
##
|
|===================================== | 52%
|
|===================================== | 53%[1] "The acceptance rate is: 0.2063"
##
|
|===================================== | 54%
|
|====================================== | 54%[1] "The acceptance rate is: 0.2062"
##
|
|====================================== | 55%[1] "The acceptance rate is: 0.2062"
##
|
|======================================= | 55%
|
|======================================= | 56%[1] "The acceptance rate is: 0.2061"
##
|
|======================================== | 56%
|
|======================================== | 57%[1] "The acceptance rate is: 0.2061"
##
|
|======================================== | 58%
|
|========================================= | 58%[1] "The acceptance rate is: 0.206"
##
|
|========================================= | 59%[1] "The acceptance rate is: 0.206"
##
|
|========================================== | 59%
|
|========================================== | 60%[1] "The acceptance rate is: 0.206"
##
|
|========================================== | 61%
|
|=========================================== | 61%[1] "The acceptance rate is: 0.2059"
##
|
|=========================================== | 62%[1] "The acceptance rate is: 0.2059"
##
|
|============================================ | 62%
|
|============================================ | 63%[1] "The acceptance rate is: 0.2059"
##
|
|============================================ | 64%
|
|============================================= | 64%[1] "The acceptance rate is: 0.2058"
##
|
|============================================= | 65%
|
|============================================== | 65%[1] "The acceptance rate is: 0.2058"
##
|
|============================================== | 66%[1] "The acceptance rate is: 0.2058"
##
|
|=============================================== | 66%
|
|=============================================== | 67%[1] "The acceptance rate is: 0.2057"
##
|
|=============================================== | 68%
|
|================================================ | 68%[1] "The acceptance rate is: 0.2057"
##
|
|================================================ | 69%[1] "The acceptance rate is: 0.2057"
##
|
|================================================= | 69%
|
|================================================= | 70%[1] "The acceptance rate is: 0.2058"
##
|
|================================================= | 71%
|
|================================================== | 71%[1] "The acceptance rate is: 0.2058"
##
|
|================================================== | 72%[1] "The acceptance rate is: 0.2058"
##
|
|=================================================== | 72%
|
|=================================================== | 73%[1] "The acceptance rate is: 0.2057"
##
|
|=================================================== | 74%
|
|==================================================== | 74%[1] "The acceptance rate is: 0.2057"
##
|
|==================================================== | 75%[1] "The acceptance rate is: 0.2057"
##
|
|===================================================== | 75%
|
|===================================================== | 76%[1] "The acceptance rate is: 0.2056"
##
|
|====================================================== | 76%
|
|====================================================== | 77%[1] "The acceptance rate is: 0.2056"
##
|
|====================================================== | 78%
|
|======================================================= | 78%[1] "The acceptance rate is: 0.2056"
##
|
|======================================================= | 79%[1] "The acceptance rate is: 0.2056"
##
|
|======================================================== | 79%
|
|======================================================== | 80%[1] "The acceptance rate is: 0.2056"
##
|
|======================================================== | 81%
|
|========================================================= | 81%[1] "The acceptance rate is: 0.2055"
##
|
|========================================================= | 82%[1] "The acceptance rate is: 0.2055"
##
|
|========================================================== | 82%
|
|========================================================== | 83%[1] "The acceptance rate is: 0.2055"
##
|
|========================================================== | 84%
|
|=========================================================== | 84%[1] "The acceptance rate is: 0.2055"
##
|
|=========================================================== | 85%
|
|============================================================ | 85%[1] "The acceptance rate is: 0.2055"
##
|
|============================================================ | 86%[1] "The acceptance rate is: 0.2056"
##
|
|============================================================= | 86%
|
|============================================================= | 87%[1] "The acceptance rate is: 0.2056"
##
|
|============================================================= | 88%
|
|============================================================== | 88%[1] "The acceptance rate is: 0.2056"
##
|
|============================================================== | 89%[1] "The acceptance rate is: 0.2056"
##
|
|=============================================================== | 89%
|
|=============================================================== | 90%[1] "The acceptance rate is: 0.2056"
##
|
|=============================================================== | 91%
|
|================================================================ | 91%[1] "The acceptance rate is: 0.2056"
##
|
|================================================================ | 92%[1] "The acceptance rate is: 0.2056"
##
|
|================================================================= | 92%
|
|================================================================= | 93%[1] "The acceptance rate is: 0.2056"
##
|
|================================================================= | 94%
|
|================================================================== | 94%[1] "The acceptance rate is: 0.2055"
##
|
|================================================================== | 95%[1] "The acceptance rate is: 0.2055"
##
|
|=================================================================== | 95%
|
|=================================================================== | 96%[1] "The acceptance rate is: 0.2055"
##
|
|==================================================================== | 96%
|
|==================================================================== | 97%[1] "The acceptance rate is: 0.2055"
##
|
|==================================================================== | 98%
|
|===================================================================== | 98%[1] "The acceptance rate is: 0.2054"
##
|
|===================================================================== | 99%[1] "The acceptance rate is: 0.2054"
##
|
|======================================================================| 99%
|
|======================================================================| 100%[1] "The acceptance rate is: 0.2054"elapsed <-proc.time() - start
print(paste0(n_it," RjMCMC iterations completed. Time elapsed:"))## [1] "1e+08 RjMCMC iterations completed. Time elapsed:"print(elapsed)## user system elapsed
## 186433.1 3.2 186463.8Process MCMC data
saveRDS(expansions, file = paste0(base_dir, "/expansions_s2.rds")) expansions <- discard_burn_in(expansions, proportion=burn_in) plot(expansions, mode="traces")
Inspect the base model parameters
We’re mostly interested in the model indicator trace and histogram, as this indicates the number of expansions. The population size trace acts as a good sanity check as well.
unique_dims <- unique(expansions$model_data$dim)
dim_counts <- sapply(unique_dims, function(i) length(which(expansions$model_data$dim==i)))
mode_dim <- unique_dims[which.max(dim_counts)]
if (mode_dim > 0) {
plot(expansions, mode="summary", k_modes=mode_dim)
} else {
plot(expansions, mode="summary")
}## Warning: `mutate_()` was deprecated in dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help## Warning: Removed 130 rows containing missing values (geom_text).
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[arrange]print(paste0("The mode for the number of expansions is: ", mode_dim))## [1] "The mode for the number of expansions is: 3"Analyse the marginal corresponding to the number of expansion mode
We will look in detail at the marginal corresponding to the mode model, and inspect the top \(K\) most likely expansions
if(mode_dim == 0) {
eval_events<-FALSE
print(paste0("No events detected"))
} else {
eval_events<-TRUE
}plot(expansions, mode="modes", k_modes=mode_dim)